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ABSTRACT 



GO 

We wish to better constrain the properties of solar flares by exploring how pa- 
rameterized models of solar flares interact with uncertainty estimation methods. We 
compare four different methods of calculating uncertainty estimates in fitting param- 
eterized models to Ramaty High Energy Solar Spectroscopic Imager (RHESSI) X-ray 
spectra, considering only statistical sources of error. Three of the four methods are 
based on estimating the scale-size of the minimum in a hypersurface formed by the 
weighted sum of the squares of the differences between the model fit and the data as 
a function of the fit parameters, and are implemented as commonly practiced. The 
fourth method is also based on the difference between the data and the model, but 
instead uses Bayesian data analysis and Markov chain Monte Carlo (MCMC) tech- 
niques to calculate an uncertainty estimate. Two flare spectra are modeled: one from 
the GOES (Geostationary Operational Environmental Satellite) XI. 3 class flare of 19 
January 2005, and the other from the X4.8 flare of 23 July 2002. We find that the 
four methods give approximately the same uncertainty estimates for the 19 January 
2005 spectral fit parameters, but lead to very different uncertainty estimates for the 23 
July 2002 spectral fit. This is because each method implements different analyses of 
the hypersurface, yielding method-dependent results that can differ greatly depending 
on the shape of the hypersurface. The hypersurface arising from the 19 January 2005 
analysis is consistent with a Normal distribution; therefore, the assumptions behind the 
three non-Bayesian uncertainty estimation methods are satisfied and similar estimates 
are found. The 23 July 2002 analysis shows that the hypersurface is not consistent with 
a Normal distribution, indicating that the assumptions behind the three non-Bayesian 
uncertainty estimation methods are not satisfied, leading to differing estimates of the 
uncertainty. We find that the shape of the hypersurface is crucial in understanding the 
output from each uncertainty estimation technique, and that a crucial factor determin- 
ing the shape of hypersurface is the location of the low-energy cutoff relative to energies 



where the thermal emission dominates. The Bayesian/MCMC approach also allows us 
to provide detailed information on probable values of the low-energy cutoff, E c , a crucial 
parameter in defining the energy content of the flare-accelerated electrons. We show 
that for the 23 July 2002 flare data, there is a 95% probability that E c lies below approx- 
imately 40 keV, and a 68% probability that it lies in the range 7-36 keV. Further, the 
low-energy cutoff is more likely to be in the range 25-35 keV than in any other 10 keV 
wide energy range. The low-energy cutoff for the 19 January 2005 flare is more tightly 
constrained to 107 ±4 keV with 68% probability. Using the Bayesian/MCMC approach, 
we also estimate for the first time probability density functions for the total number 
of flare accelerated electrons and the energy they carry for each flare studied. For the 
23 July 2002 event, these probability density functions are asymmetric with long tails 
orders of magnitude higher than the most probable value, caused by the poorly con- 
strained value of the low-energy cutoff. The most probable electron power is estimated 
at 10 28,1 erg sec -1 , with a 68% credible interval estimated at io 28 - 1_29 -°erg sec -1 , and a 
95% credible interval estimated at io 280-30 - 2 erg sec -1 . For the 19 January 2005 flare 
spectrum, the probability density functions for the total number of flare accelerated 
electrons and their energy are much more symmetric and narrow: the most probable 
electron power is estimated at io 27 - 66±0 - 01 erg sec -1 (68% credible intervals). However 
in this case the uncertainty due to systematic sources of error is estimated to dominate 
the uncertainty due to statistical sources of error. 



Subject headings: Sun: flares 
- methods: statistical 



Sun: X-rays, gamma rays — methods: data analysis 



1. Introduction 



The detailed understanding of solar flares requires an understanding of the physics of ac- 
celerated electrons, since electrons carry a la r ge fra ction of the total energy released in a flare 
(|Lin &: Hudsonl Il97ll . 119761 ; lEmslie et al.l [2004J, 120051 ) . Since we cannot measure the electron flux 
in situ, the behavior of the flare-accelerated electrons is inferred from the photons emitted by 
their interaction with the ambient plasma. For a general inhomogeneous optically thin source of 
plasma density n(r) and electron flux densitjo energy spectrum F(E, r) (electrons cm -2 s -1 keV - ) 
in volume V for electron energy E, the bremsstrahlung photo n flux density energy spectrum 1(e) 



(photons cm 2 s 1 keV at Earth distance R) can be written ( Brownlll97ll : lBrown et al.ll2003l ) as 
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nV 
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F{E)Q{e,E)dE, 



(1) 



x In this paper, "flux density" refers to an amount per unit area per unit time. 



where n = J v ndV/V, F(E) is the mean electron flux distribution, F(E) = J v n(r)F(E, r)dV/(nV), 
and Q(e,E) is the bremsstrahlung cross-section differential in photon energy e. In this paper we 
model the photon flux density energy spectrum as the sum of emission due to a flare-injected 
electron flux spectrum interacting with a target, and emission from hot plasma with a Maxwellian 
distribution of speeds corresponding to some temperature T. 

The Ramaty High Energy Solar Spectroscopic Imager (RHESSI, Lin et al. 2002) flags all 
photons detected in any one of the nine germanium detectors by the time of occurrence (to 1 
microsecond), the amount of energy lost by the photon in the detector (in 0.3-keV-wide pulse 
height analyzer (PHA) bins), and the detector number. For spatially integrated spectral analysis, 
the counts can be combined arbitrarily over different detectors and PHA bins. 

We define D= (-Di, ..., Di, ...D n ) as the number of counts observed in a given set of energy-loss 
bins labeled in the range 1 < i < n in a given time interval. These counts are noisy, and are 
assumed to be drawn from a Poisson distribution with a mean of Cj, 

P ( Dl ) = ^-e- c >. (2) 

The measured count rate Rf in energy-loss bin T is determined from the measured counts Di 
divided by the live timqj t^x- The predicted count rate Rf arises from the incident photon flux 
rate via 

R? = M l3 I r , (3) 

that is, the predicted count rate in an energy-loss bin T is modeled via a detector response matrix 
Mij for an incident photon flux spectrum Ij, where the index j, 1 < j < m, labels energies at which 
the incident photon spectrum is calculated. The response matrix M^ is calculated by RHESSI 
Solarsoft routines once the count energy-loss bins (indexed by 'i') and incident photon energies 
(indexed by 'j') are defined. The incident photon flux energy spectrum is deduced by comparing 
the observed with the predicted count rates in all energy bins assuming a model for the photon 
flux energy spectrum until some criterion for agreement is met. 

One goal of RHESSI data analysis is to recover the electron flux energy spectrum F(E, r) from 
the detected counts Di in a given time interval. In general, this requires detailed knowledge of 
the energy losses suffered by the bremsstrahlung-producing electrons in the emitting volume. It is 
often only practical to recover F(E); to do this, two approaches are commonly taken. 

Since the rates are measured, and everything other than F(E) is known (either calculated, 
measured or assumed), F(E) can be obtained through Equations Q] and [3l This approach is known 
as inversion. The advantage of inversion is that one does not make an assumption as to the nature 
of the mean electron flux distribution. The disadvantage of this approach is that noise in the 



2 The live time is the observation time minus the dead time. The dead time is the amount of time that the detector 
cannot respond to an incoming photon. 



observed data and errors in instrument calibration can lead to the creation of spurious features in 
the solution. This effect can be mitigated by adding extra constraints to the inversion process which 
forces the solution to be smooth across energy bins (note that this is required by the bremsstrahlung 
process and RHESSI's energy resolution). Consider discretizing Equation Q] by energy bins to yield 
a matrix expression, 

J = AF (4) 

where J is a m-element vector representing the observed number of photons /(e), A is a m x 
n=matrix representing Q(e, E) and F is a n-element vector representing the mean electron spectrum 
F(E). The standard approach is to minimize the residual 

IIAF- Jll 



for F where || • || is the Euclidean norm. This matrix problem can be ill-posed due to the noise 
sources discussed above, or by A being ill-conditioned or singular. Regularization mitigates these 
issues by imposing extra constraints on the solution for F. Tikhonov regularization does this by 
adding an extra term ||rF|| for some choice of Tikhonov matrix T, to the above minimization 
problem, yielding 
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IAF- Jll + ||rF| 



(5) 



Piana et al.1 (|2003l ) demonstrate a Tikhonov-regularized inversion algorithm that takes the observed 
counts and finds F(E) and the uncertainty on F(E). iPiana et al.l (120031 ) show that this method 
led to an unexpected 'dip' in the mean electron spectrum which is thought (in most cases) to arise 
from the presence of a sig nificant photospheric albedo flux contributing to the observed X-ray flux 
(JKontar et alJbood . bpO^ . 



In the second approach, known as forward fitting, a parameterized model for the mean electron 
flux distribution F(E) is used to describe the photon flux Ij incident at RHESSI. The photon 
emission, parameterized by 6 (Ng variables) is 



Ij(8)- 



(6) 



A fitting process is then used to find values of the parameters that best reproduce the counts Di 
observed by RHESSI. The disadvantage of this method is that the spectral model is prescribed 
rather than derived, and so features that are not in the mode l cannot be describ ed by it, although 
their presence in the data may be indicated by the residuals (JBrown et al.ll2006l ). The advantages 
of this method are that by judicious choice of parameterization the major features of the spectrum 
can be modeled, and values to the parameters with uncertainty estimates can be obtained. 

In this paper, we use the forward fitting approach and consider four different methods of 
estimating a range of "acceptable" model parameter values that describe our understanding of the 
flare within the confines of the model. By comparing different methods, we seek to understand the 
differences in the final answer that may be brought about by the way the estimates were obtained. 
Further, by comparing two different spectra we can better understand how, for a given model, the 



estimated parameter values and errors are influenced by the data. It is assumed that the only 
source of noise is the Poisson distribution that follows naturally from independent photon events 
(Eq. ED. 

Systematic err or sources are u ndoubtedly important in determining the uncertainties in the 
model parameters ()Lee et al.ll201ll ). but they are not explicitly included in the uncertainty deter- 
mination methods described below. Two types of systematic uncertainties are common in this type 
of spectroscopy, integral and differential. Integral uncertainties are basically the uncertainties in 
the overall sensitivity of a given detector. Based on comparisons of flare spectra measured with 
different detectors, they are known to be smaller than approximately 10%. They affect primarily 
the absolute value of the emission measure in the thermal model and the total electron flux in the 
nonthermal electron spectrum. The differential uncertainties are basically the uncertainties in the 
sensitivity in each energy bin with respect to its neighbors. They affect primarily model parameters 
that depend on the slope of the measured spectrum. They are therefore important for the tem- 
perature in the thermal model and the low energy cutoff and power-law index of the nonthermal 
electron spectrum. For RHESSI, the differential un certainties are less than 1 % and are generally 
negligible compared to the statis t ical u ncertainties. iMilligan h Dennis! (J2009I ) (using detectors 1, 
3, 4, 5, 6, and 9) and ISu et al.l (|201ll ) (using detectors 1, 3, 4, 6, 8 and 9) show that there is 
scatter in the best-fit parameter values determined from different individual detectors for the flare 
models they considered but that the range of the scatter indicates that the systematic errors are 
not significantly greater than the statistical errors. The systematic uncertainties are not important 
in developing a basic understanding of how each uncertainty determination method behaves in the 
presence of noisy data and consequently they have not been included in the analysis done for this 
paper. 

The observations and spectral models are described in Section [2l Section [3] describes the 
parameter and uncertainty estimation methods used. Section U] describes the results and Section [5] 
discusses the implications of these results for fitting spectral models to RHESSI data. 



2. Spectral model and observations 



In the X-ray energy range covered by RHESSI (JLin et al.ll2002l ) - generally from ~3 keV up to 
a few hundred keV - the emitted photon spectrum is modeled as the sum of a thermal component 
that generally dominates at the lower X-ray energies, typically below ~10-20 keV, and a non- 
thermal component that dominates at higher energies. The thermal component is the line and 
continuum emission from the flare-heated plasma. The line emission is mainly from transitions in 
highly ionized iron - primarily FeXXV - that appears in the RHESSI spectrum as an unresolved 
peak at 6.7 keV with a much weaker feature at ~8 keV. The continuum emission is a combination 
of free- free emission (bremsstrahlung) and free-bound emission (recombination radiation). 



For our spectral analysis, we have used the thermal line-plus-continuum spectra provided by 



C HIANTI (IDere et al.lll997l . l2009l ) assuming an isothermal plasma with the ionization balance given 
bv lMazzotta et al.l (|1998l ) and the "sun coronal" abundances given bv lFeldman et al.l (J1992I ). The 
only free parameters are the temperature (kT in keV) and the volume emission measure (EM in 
cm" 3). 



The thermal continuum emission is made up of the sum of bremsstrahlung (or free-free) emis- 
sion and free-bound emission. The form of the bremsstrahlung contribution as a function of photon 

energy e is approximately 

\rr,M] 

-e/kT), (7) 
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wher e k is Boltzmann's constant and Ithermai is in units of photons s *erg 1 ( Tandberg-Hanssen &; Emslie 
1988I ). The free-bound continuum spectrum has a similar dependency on EM and T. 



The non-thermal component of the measured X-ray spectrum is bremsstrahlu ng from flare- 
accelerated electrons interacting with the ambient medium. Following iBrownl (|197ll ) , we assume a 
cold, thick target, meaning that the electrons collisionally lose their energy in cold, fully ionized 
plasma as they radiate. The energy loss rate per unit distance x as an electron with speed v streams 
through the ambient plasma is dE/dx = —2Kn e (x)/(mv 2 ), where m is the el ectron mass, n e (x) is 
the number density of plasma electrons, and K is approximately constant (see lHolman et al.ll201ll ). 
Using this result, the mean electron flux becomes 



F (E) = - 



1 mv 2 



nV 2K 



F (E )dE , 



where Fq(Eq) is now the injected electron flux energy spectrum (electrons s l keV 
following broken power-law functional form for the spectrum of injected electrons: 



(8) 



We use the 



Fo(Eq) = A < 





(E /E p )^ 

(E /E p y s -(E h /E p ) 5 ^ 





Eq < E c 
E c < Eq < Ef, 
Eb < E < Eh 
Eq > Eh 



(9) 



The seven parameters of this nonthermal component are the normalization parameter A, the low- 
and high-energy cutoffs, E c and Eh, the pivot energy E p , the break energy Eb, and the power-law 
indices below and above the break energy, 5\ and 82 , respectively. The radiated X-ray spectrum is 
modeled as the sum of the isothermal component and Equation[H where F(E) is given by Equations 
[8]and[9j The X-ray emission is assumed to be isotropic and, with this assumption, the contribution 
flux from phqtosph eric albedo to the total incident X-ray at the instrument can be estimated (see 



Kontar et al. 



2011 



We model RHESSI spectral data from two flares - the GOES class XI. 3 flare on 19 January 
2005 starting at 08:03 UT, and the X4.8 flare starting at 00:18 UT on 23 July 2002. We choose 
these flares because previous studies have shown that the low-energy cutoff - E c - is estimated to lie 



in very different portions of the spectrum. In the 23 July 2002 event, the low-energy cutoff of the 
flare-accelerated electrons is estimated to have an energy in the region where the observed hard X- 
ray emission is thermally dominated. This makes it difficult to place limits on the low-energy cutoff 
since it is difficult to determine the signal of the flare-accelerated electrons against the dominant 
thermal bremsstrahlung emission. Most flares are thought to have low-energy cuto ffs close to or in 



the re gion where the emission is dominated by thermal bremsstrahlung. In contrast. I Warmuth et al. 



(J2009I ) studied the 19 January 2005 event, and found that late in the impulsive phase, the low- 
energy cutoff energy much higher than energies at which the thermal bremsstrahlung dominates. 
Therefore, thermal bremsstrahlung cannot be a significant factor in determining the uncertainty in 
the low-energy cutoff for this flare. The low-energy cutoff is one of the most important properties 
of a flare as its value strongly influences the estimated flare-accelerated electron energy content. 
Therefore, knowledge of the uncertainty in the low-energy cutoff directly influences knowledge of 
the energy content of the flare. Hence, these two flares and the models used to study them are good 
test-beds for understanding how different uncertainty estimation methods operate when generating 
uncertainties for parameters that are crucial for understanding the properties of solar flares. 

Table [1] has details of the two flares and the two spectral accumulation times chosen, the 
models used, and the best-fit parameter values obtained that fit the spectral models to the data 
(see Section I3.ip . These two spectra were chosen because they were both well observed with 
RHESSI and they highlight the excell ent spectral capabilities of the cooled germanium detectors 
of this instrum ent (ISmith et al.l 2002 ). Both flares have been extensive ly analyzed previous ly - 
see for example I Warmuth et al.l (120091 ) for the 19 January 2005 flare and lHolman et al.l (120031 ) for 
the 23 July 2002 flare. The most notable difference between the two spectra is that the first has 
a low-energy cutoff in the electron spectrum of over 100 keV, well above the thermal component. 
This is in contrast to the s econd flare where the low-energy cutoff is estimated to be below ~ 
40 keV (JHolman et al.ll2003l ) and consequently difficult to determine because of the dominance of 
the thermal component at lower energies. This difference between these two flares motivates their 
selection for this study. These two flare events are good candidates that allow us to explore how 
well we can determine the value of the crucial low-energy cutoff parameter (and flare properties 
that depend on it) given the data, the model, and the uncertainty estimation methods used. 

Traditionally, RHESSI spectral analysis invol ves summ i ng da ta from multiple RHESSI detec- 
tors to improve counting statistics - see for example lSu et al.l (120091 ). Instead of this usual approach, 
we chose to use data from just one detector with good energy resolution and sensitivity - detector 
#4. This allowed us to apply the most accurate corrections for energy resolution and calibration, 
pulse pile-up, and background subtraction. In the time periods selected, the count rates were suf- 
ficiently high that selecting a single detector did not seriously degrade the spectroscopy capability 
up to the highest energies considered of ~500 keV. The energy bin widths were chosen to be as 
narrow as possible to preserve spectral details resolvable with the detector's ~1 keV FWHM spec- 
tral resolution while maint aining >30 c ounts in each bin as required for the x 2 analysis procedure 
to be approximately valid IWassermanl ()2003l ). The only part of the spectral data that is affected 



by small numbers is at the high energy part of the spectrum, well away from the low energy part 
of the spectrum. At these energies, the simple Normal approximation to the Poisson distribution 
- (Poisson(X) m N(\, y/\) for A 'large') - is no longer appropriate. However, the gross properties 
we are most interested in - flare energy content, the number of flare-accelerated electrons and the 
probability density function of the low-energy cutoff, are largely unaffected by a biassed fit of a 
spectral model to the data at high energies, since these properties are largely determined by the 
flare spectrum at energies where the Normal distribution can be used. We can assert this for the 
flares studied in this analysis because these are relatively large flares with large numbers of counts. 
The vast majority of flares are smaller than the ones studied here, and therefore fits or parame- 
terized models to the data are more likely to suffer from biassed fits over more extensive energy 



ranged 3 ' 



Both flares have been extensive ly analyzed pr e yiousl y - see for example IWarmuth et al.1 ()2009l ) 



for the 19 January 2005 flare and iHolman et al.l (J2003I ) for the 23 July 2002 flare. For ease in 



comparing results in each case, we have generally followed their lead in choosing background spectra, 
energy ranges, model components, fitting procedures, etc. in the spectral analysis. Table Q] has 
details of the two flares and the two spectral accumulation times chosen, the models used, and 
the best-fit parameter values obtained that fit the spectral models to the data (see Section [3~T]) . 
Corresponding count fluxj and photon flux spectra are shown in Figures Q] and [2j The model 
count flux spectrum is computed by taking the best fit photon spectrum and convolving it with 
the instrument response matrix. Figures [lb and [2b shows the best fit photon spectrum and the 
photon spectrum derived from the measured count flux spectrum using the ratio of the best fit 
photon spectrum to the measured counts in each energy bin. The units in Figures [T] and O) are 
photons s^ 1 cm" 2 keV^ 1 .) 



2.1. 19 January 2005 

The first flare considered was the GOES XI. 3 flare that peaked at 08:22 UT on 19 January 
2005 on the solar disc at N15W51. We used the RHESSI observations of this flare from 08:26:00 - 



08:26:20 UT, the same time interval when lWarmuth et al.l (|2009i ) found an unusually hard spectrum 



during the final peak of the impulsive phase, possibly resulting from a low-energy cutoff in the 
electron spectrum as high as 120 keV (see their Figure 1 for RHESSI light-curves of this event). 

We used the standard procedures that form OSPEX, the standard spectral analysis package 
used in RHESSI data analysis, to determine the best-fit parameters of the thermal and nonthermal 



3 It should also be noted that even although a substantial part of the spec trum have la rge enoug h counts, biassed 
values to the fit are still possible when minimizing a x 2 -hke expression - see ICashl ( f 979 ) and also iHumphrev et al 
1 2009T ) and references therein. 



4 By count flux we mean the measured count rate per keV divided by a nominal detector area corrected for grid 
transmission, equal to 38cm 2 for the single detector used in our analysis. 



components of the incident photon spectrum. As is common in RHESSI data analysis, the back- 
ground spectrum that was subtracted from the measured count rate spectrum was calculated by 
linear interpolation in time between spectra measured before and after the flare. The estimated 
background spectrum is about an order of magnitude less than the flare spectrum at all energies 
considered. The background can therefore be considered as having very little influence on the 
final probability density functions of the model parameters and the gross proper ties of the flare 



such a s its energy content and the number of flare accelerated electrons. Following IWarmuth et al. 



(|2009l ). we included two narrow Gaussian-shaped emission lines in the model photon spectrum to 



accommodate features in the count-rate spectra that are believed to be instrumental in origin. 

We included the standard corrections for energy calibration adjustments and pulse pile-up, 
but these did not play a significant role for the selected time interval since the attenuators were in 
the A3 state (both thick and thin attenuators in place) resulting in relati vely low counting rates . 



The albedo component was not included here, although it was included bv lWarmuth et al,l (|2009l ). 
We found that adding the albedo component did not significantly alter the fitted parameters or 
the estimates of the uncertainties. We used the following energy bins for this flare: 1/3 keV from 
3 to 15 keV, 1 keV from 15 to 50 keV, 5 keV from 50 to 100 keV, and 10 keV from 100 to 300 
keV. The photon spectrum was extended above the fitted energy range up to 600 keV to allow for 
non-photopeak response of the detector. 



Again, following IWarmuth et al.l (120091 ). we modeled the therma l component wi t h a s ingle 



temperature function from CHIANTI using coronal abundances and a lMazzotta et al.1 (119981 ) ion 



ization balance. The nonthermal component was modeled assuming thick-target interactions of 
electron with a single-power-law spectrum at energies above E c . This is accommodated in Eq. [9] 
by fixing both 62 at the default value of 6.0 and Et, at 32 MeV so that they have no significant 
effect on the bremsstrahlung X-ray spectrum in the fitted photon energy range below 300 keV. E^ 
was fixed at 32 MeV so that, like Ef,, it has negligible effect on the bremsstrahlung X-ray spec- 
trum in the fitted energy range, and so is equivalent to having no cutoff at all. For this flare, the 
normalization was taken to be Fq, the total integrated electron flux over the electron energy range 
from E c to Eh with E p fixed at 1 keV, instead of A in Eq. [9l The advantage in normalizing to Fq 
is that this is a physically interesting quantity. The disadvantage is that it is strongly dependent 
on the value of both the low-energy cutoff and the spectral index. For the conditions described 
here, Fq = AE^~ /{5\ — 1). The package OSPEX was configured to use this implementation of 
Equation [9] for this flare. An alternate implementation was required for the 23 July 2002 event (see 
Sections I2T21 and IO). 



In our detailed spectral analysis and assessment of uncertainties, we had a total of seven free 
parameters - EM, kT, Fq, E c , Si, G\ and G2 ~ (see Table dj). Other parameters covering the 
instrumental effects - energy calibration, pulse pile-up, and Gaussian features below 10 keV - were 
determined from the analysis of the count-flux spectra for other time intervals and other flares, and 
then fixed for the subsequent determination of uncertainties in this time interval. The amplitudes 
of the two Gaussians (G\ and G2) were free during the spectral fits. 
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(a) Count Flux Spectrum, 19-January-2005 



(b) Photon Flux Spectrum, 19-January-2005 
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Fig. 1. — Count and photon spectra for the 19 January 2005 flare in the analysis period 08:26:00 
to 08:26:20 UT. (a) The histograms with ±1<t statistical error bars represent the background- 
subtracted count fluxes (black) and the background fluxes (pink) vs. energy loss in the detector. 
The smooth curves represent the different components of the model used to fit the data as follows: 
isothermal (green), thick-target bremsstrahlung (yellow), Gaussians (blue and cyan). The sum of 
all the components is shown in red. (b) Incident photon flux (in units of photons s cm~ 2 keV ) 
vs. photon energy with the different components of the model shown in the same colors as in (a). 
The energy range used for the spectral fits lies between the vertical line at 6.45 keV and the edge 
of the plot at 300 keV. 
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2.2. 23 July 2002 



The second flare considered was the GOES X4.3 flarqfl that peaked at 00:35 UT on 23 July 



2002 from a location closer to the limb at S13E72 than the first event. Following iHolman et al. 



(J2003I ) . we chose to analyze the time interval from 00:30:00 to 00:30:20.250 UT during the first 



pea k of the impulsiv e phase (see their Figure 1 for RHESSI X-ray light-curves of this event; see 



also iLin et al.1 (120031 ). their Figure 1 for a lightcurve of the GOES X-ray flux). The measured 
X-ray spectrum was again assumed to be the sum of an isothermal spectrum and the thick-target 
bremsstrahlung spectrum from non-thermal electrons with the broken power-law of Eq. EO In this 
case, the full double power-law was assumed with the break energy, Eb, and the second power-law 
index, 62, both free parameters. The normalization constant for this flare, A in Eq. [9J was defined 
as the electron flux at the pivot energy E p that was fixed at 50 keV. As with the first flare, the high 
energy cutoff to the electron spectrum Eh was set at 32 MeV to ensure that it had no significant 
effect in the fitted photon energy range. 

The following 130 energy bins were used for this event: 1-keV wide bins from 3.0 to 40 keV, 
3-keV from 40 to 100 keV, 5-keV bins from 100 to 150 keV, 10-keV bins from 150 to 500 keV, 
1-keV bins from 501 to 520 keV, and 10-keV bins from 520 to 600 keV. We extended the energy 
range of the assumed photon spectrum up to 20 MeV to allow for the off-diagonal elements of the 
instrument response matrix due to the non-photopeak response of the detector. The fitted photon 
energy range was restricted to be above 15 keV to avoid the need for the two Gaussian emission 
line sources to accommodate the supposed instrumental features below 10 keV used for the first 
flare. The upper energy of the fit range was extended up to 500 keV to provide more information 
on the power-law spectrum above the break energy. This increase in the upper energy limit also 
necessitated a dding in a nuclear co mponent in the form of a template appropriate for a power-law 



ion spectrum (jMurphy et al.lll991l ) with the normalization parameter fixed at the value obtained 
to give a best fit to the data. This nuclear component (shown in Fig. [2]) contributes <10% to 
the photon flux at all energies below ~400 keV and hence has only marginal significance in the 
subsequent analysis. 

Other parameters were determined from least-squares fits to the count-flux spectrum and then 
fixed for the subsequent determination of uncertainties. These included parameters to characterize 
the instrumental effects of pulse pile-up that is a more important component for this flare since the 
count rates were a factor of ~10 higher than in the first flare. Also, although it is not significant for 
flares at the solar limb, the albedo spectrum was included f o r this flare assuming isotropic X-ray 



emission using the the procedure described in iKontar et al.l (J2006I ) and implemented in OSPEX. 



Both the pile-up and albedo components are shown in Fig. [2l 

For our detailed spectral analysis and assessment of uncertainties for this flare, there was a 



5 Many more details concerning this flare can be found in the special issue of the Astrophysical Journal Letters 
(vol. 595) dedicated to its study. 
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total of seven free parameters - EM, kT, A, E c , E^, S\, 82 (see Table[]jJ. The background-subtracted 
count flux and photon spectra are shown in Fig. [2] along with the best-fit model components. Note 
that the implementation of the normalization used for this analysis is different from that used for 
the 19 January 2005 flare. In this analysis using the normalization A at the pivot energy E p is 
preferred. The reason for this choice is given in Section [4.2i 



a) Count Flux Spectrum, 23-July-2002 



1.0000=- 



0.1000 7 



0.0100 7 




(b) Photon Flux Spectrum, 23-July-2002 




Energy (keV) 



29-Nov-2011 18:52 



Energy (keV) 



30-NOV-2011 17:21 



Fig. 2. — Similar to Fig. Q]for the 23 July 2002 flare. The following three additional components are 
included in this plot: albedo (purple), pulse pile- up (blue), and the nuclear template (cyan). The 
two Gaussians shown in Fig. Q]were not used for this fit. The energy range used for the spectral 
fits lies between the two vertical lines at 15 and 500 keV. 



3. Parameter and Uncertainty Estimation Methods 



Four different methods of uncertainty estimation are described below. The first three methods 
- 'covariance matrix', 'x 2 -mapping' and 'Monte Carlo' sampling (Sections 13.1. 1| 13.1.21 and 13.1.31 
respectively) are widely used to estimate errors in parameter values. The fourth method is based 
on Bayesian probability and the Markov chain Monte Carlo (MCMC) method fSection l3.2.ip . Each 
of these methods is applied to the spectral model and data described in Section [5J and the results 
are tabulated in Table [2] (19 January 2005) and Tabled (23 July 2002). 
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Table 1. Flare characteristics and model parameters. 









Flare 1 




Flare 2 


Date 




19 January 


2005 


23 July 2002 


GOES Start/Peak/End Times 


08:03/08:22/08:40 UT 


00:18/00:35/00:47 UT 


GOES Class 






X1.3 






X4.8 


Location on the Sun 




N15W51 




S13E72 


Radial distance from Sun center 1 




763" 






904" 


Time Interval Analyzed 




08:26:00 


- 08:26:20 UT 


00:30:00 


- 00:30:20.250 UT 


Fitted Photon Energy Range 


6.45 to 300 keV 


15 to 500 keV 


Fitted Photon Energy Bins 






90 






90 


Parameter 


Units 


Value 2 


Free/Fixed 3 


Value 2 6 


Free/Fixed 3 


Thermal Plasma 














EM 


10 49 cm" 3 


2.31 




free 


2.16 


free 


Temp. (kT) 


keV 


2.03 




free 


3.18 


free 


Abundance 


coronal 


1 




fixed 


1 


fixed 


Non-thermal Electrons 














Fo, integrated flux 4 


10 35 s" 1 


0.17 




free 




not used 


A, flux 5 at E p 


10 35 s^ 1 keV -1 


not used 


0.028 


free 


E c 


keV 


105 




free 


32.0 


free 


E p 


keV 


1 




fixed 


50 


fixed 


E h 


keV 


32,000 




fixed 


256 


free 


E h 


keV 


32,000 




fixed 


32,000 


fixed 


Si 




3.57 




free 


3.40 


free 


s 2 




6.0 




fixed 


3.92 


free 


Nuclear Template 














Normalization 


photons cm -2 


not used 


2.11 


fixed 


Gaussians 














Gi peak E 


keV 


8.44 




fixed 




not used 


G2 peak E 


keV 


9.95 




fixed 




not used 


Gi amplitude 


photons cm- 2 s" 1 


33,300 




free 




not used 


Gi amplitude 


photons cm- 2 s -1 


12,800 




free 




not used 


Gi, 2 FWHM 


keV 


0.1 




fixed 




not used 



x As measured in the heliocentric-cartesian (heliographic) co-ordinate system ( Thompsonll2006l ). 

2 Best-fit value of parameter computed using OSPEX - see Section |3. II 

3 Parameter fixed or allowed to go free in OSPEX least-squares fitting. Parameters noted as 'fixed' are frozen at 
their values in subsequent uncertainty analyses. 

4 Total electron flux integrated over all energies from E c to Eh- 

5 Electron flux at E p . 

6 The use of the pivot value in the implementation of Equation [5] is explained in Sections 12. II and 12.21 
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3.1. Methods 1-3: Parameter and uncertainty estimation via nonlinear 

least-squares fitting 

The first three methods are based on finding a local minimum Xmin ^° the quantity 

X 2 = E ~ 2 ] • (10) 

1=1 ' 

for some value of = and i«j. The quantity x 2 is a hypersurface parameterized by 0. The 
quantity # is found by performing a nonlinear weighted least squares fit minimizing x 2 with re- 
spect to 0. There are many different ways of implementing this minimization. The minimization 
was achieved using the OSPEX spectral analysis package which uses the IDL/Solarsoft routine 
MCURVEFIT.pr o. This routine is based on the nonlinear least-squares Levenburg-Marquardt fit- 



ting algorithm of lPress et al.l (J1992I ) (pages 675-683). This implementation of the algorithm ignores 



the second derivative of the fitting function R i with respect to 0, and is therefore equivalent to 
assuming that the fitting function is linear with respect to near the best-fit value 0. 

The value of is derived as follows. The process is begun with an initial estimate of 0, 0°. The 
corresponding flux rate spectrum R { is calculated and W{ is set to \/Ci{0 G )/tLT- This value of 
u>i is passed to MCURVEFIT.pro. This routine refines the estimate of the values of the spectral 
parameters, stopping when the termination condition is mefl This first estimate is to is labeled 
1 . The fitting routine is run again this time using 1 as the initial estimate to and with wi set 
to y/Ci(W)/tLT, yielding a second estimate 2 . The routine is run a third and final time using 2 
as the initial estimate to and with u>i set to \/Ci(0 2 )/t[ J T, yielding a final parameter estimate, 
labeled 0. 

Estimates of the uncertainty in the value are found by defining a scale-size of variation in 
the x 2 -hypersurface around in different ways. Three different methods of defining and estimating 
the uncertainty in the value are described below. 



3.1.1. Method 1: Uncertainty Estimation by Estimating the Covariance Matrix 

This method uses the curvature matrix of the v^hypersurface evaluated at to estimate the 
uncertainty in each parameter, via the assumptions that the measurement errors in the data D are 
Normally distributed and that either the model R- is linear in its parameters, or the region over 
which the uncertainty estimate spans can be replaced by a linear approximation to the original 
model. 

The curvature matrix a of the the x 2 -hypersurface arises in linear and nonlinear least-squares 



6 MCURVEFIT.pro stops iterating the Levenburg-Marquardt fitting algorithm when the relative change of x 2 from 
its current value to its previous value is less than 0.001. 
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fitting algorithms and is denned as a^ = d 2 (x 2 )/(d9id9j) for 1 < i,j < Ng. The implementation 
of MCURVEFIT. pro gives an unce rtainty estimate to each of the free parameters based on the 
curvature matrix (jPress et al.lll992l ). The uncertainty for 9i (for 1 < i < Ng) is 



dOi 



± 



(11) 



when evaluated at 9 = 9 (the value that minimizes x 2 > Equation I10p . The quantity a -1 in the 
right-hand side of Equation [11] is the matrix inversion of the curvature matrix and is an estimate of 
the covariance matrix of the fit parameters, evaluated at 6. Its diagonal elements are the covariance 
scale-sizes that defines th e uncertaint y estim ates in this method. Full details of the derivation of 
Equation Qj] are given in iPress et al.l ()1992l ). pages 690-692. The assumptions in this derivation 
also imply that the probability distribution for 80 o b s (the expected error in the value of 0) is a 
multivariate Normal distribution around 9. The uncertainty estimate given by Equation QT] is 
quoted as the 68% value in Tables [2] and [3] 



3.1.2. Method 2: Uncertainty Estimation using x 2 -mapping 



In this method, parts of the shape of the x 2 -hypersurface around Xmin are explicitly calculated. 
It is assumed that the value of the x 2 -hypersurface as defined by Equation 1101 at a particular point 
9, follows a ^-distribution. By fixing a probability and finding where that probability occurs as 
a function of the parameters, one can measure scale-sizes in the x 2 -hypersurface that define an 
estimate of the uncertainty in the value of 9 with that probability. The procedure is described 
below. 

One of the parameters 9 in the set 9 is stepped through a range of values while the others are 
allowed to vary so as to minimize x 2 , yielding a value y 2 • The quantity 5y 2 = x\ ~Xmin * s assumed 



to have a ^-distribution with one degree of freedom ( Press et al.lll992l ). For such a distribution 



one can therefore expect that 5x 2 < 1 occurs approximately 68% of the time and 5x 2 < 4 occurs 
approximately 95% of the time. Values for the 68% and 95% confidence intervals are found where 



tx 2 



3 68% 



95%- 



) = l,5 X \9^) 



(12) 



respectively. The uncertainty estimates defined by this method are quoted as differences from 9 in 
Tables M and El that is, 

9]™< 1% -9 i (13) 

for 1 < i < Ng where q = 0.68 and q = 0.95 and 9 i q ° is defined by Equation [12j Typically there 
are two values of 9 i q ° that satisfy Equation [12] corresponding to the upper and lower confidence 
limits of the parameter value 9{. When no value of 9 can be found that satisfies the conditions of 
Equation I12| this is reported as 'not determined in Ta bles [2] and [3] Fina lly, this method uses the 
same underlying assumptions as those in Section [3.1.11 (JPress et al.lll992l ). 
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3.1.3. Method 3: Uncertainty Estimation using the Monte Carlo method 

This method of obtaining uncertainty estimates on 9 is commonly called the "Monte Carlo" 
method. This method begins by assuming that the value 9 found in method 1 best describes the 
observation via the parameterized model. By Equation [3l this defines an estimated count flux rate 
spectrum of R, that is assumed to be a good estimate of the true count flux spectrum. Estimates 
of the errors in 9 are found by generating a new spectrum such that counts in energy-loss bin i are 
drawn from a Poisson distribution with mean value i?^ for all 1 < i < n. This new spectrum is 
now fit using the same physical model and fit process as the original fit generating 9. The sampling 
and fitting process is repeated; the distribution of values found is centered at 9, and the width of 
distribution estimates the uncertainty in 9. The sample and fit process is repeated 10,000 times, 
from which normalized frequency distributions F(9i) (1 < i < Ng) are calculated. The uncertainty 
estimate used excludes the tail values in a frequency distribution F(9i). The W0q% uncertainty 
estimate for < q < 1 is defined as [# L |q,#^|q] where 

)L \ q r°° 



F(9i)d9i = / F(6i)d9i = (1 - q)/2. (14) 

Je H \ q 

This definition finds an interval [9 \ q ,9 \ q ] such that 100(/% of the measurements are within the 
interval and an equal percentage of the measurements are both above and below the interval. This 
definition of the interval is also guaranteed to contain the median value (which can be found from 
Eq. HH by setting q = 0) . The uncertainty estimates found by this method are quoted in Tables [2] 
and [3] as differences 



^ L \ q -9 l ,e H \ q -9 l (15) 



for q = 0.68 and q = 0.95 (1 < % < Ng). 



3.2. Method 4: Parameter and uncertainty estimation using Bayesian data analysis 

This method uses parameter and uncertainty estimation based on Bayesian data analysis meth- 



ods (I Jaynesl 120031 : lGregoryll2005l ). In Bayesian data analysis, the probability of a hypothesis H is 
calculated via Bayes' theorem. Denoting by p(a\b, c) the conditional probability that proposition a 
is true given that propositions b and c are true, Bayes' theorem is 

p( „ |DiI) ^jwp) (16) 

where H is the hypothesis to be tested, D is the observation, and I is any other applicable infor- 
mation we have prior to calculating the posterior. 

The left hand side p{H |D,X) is called the posterior probability of the hypothesis, given the 
data and the prior information, and it encapsulates the available knowledge about the hypothesis. 
The quantity p{H\X) is called the prior distribution and represents what we know about H prior to 
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calculating the posterior. Often a prior describes a probability density function of likely parameter 
values. The sampling distribution or likelihood, p(D\H,I), represents the likelihood of the data 
given the hypothesis H and information X. The quantity p(D\I) is the unconditional distribution 
of D and is a constant which ensures that the posterior integrates to 1. 

In this paper, the hypothesis H is that a model count spectrum C parameterized by 9 explain 
the observations D. Since the counts in each energy bin are Poisson distributed, the likelihood of 
measuring a certain set of counts Ci(9 B ) becomes 

p(n\9 B ,l) = H^y-e- c ^ B l (17) 

t=i *' 

Each parameter in the fit has its own prior p(0j.|Z), 1 < k < Ng so that p{H\X) = JJfLi P(@k \Z) ■ 
Each parameter is given a flat or uniform prior in a fixed range, that is, there is an equal probability 
that the parameter can take any value in the fixed range. Table d] tabulates the permitted range 
of values for each parameter for each model. 

The Bayesian posterior probability that a set of values 9 explains the observations D is 
proportional to the product of the likelihood and the prior. The posterior summarizes the complete 
state of knowledge of 9. Values that give rise to higher posterior probability are better explanations 
of the data, and vice versa. The best explanation of the data is the maximum a posteriori (MAP) 
value 9 M which maximizes the value of the posterior. Under the Bayesian interpretation of 
probability, values 9 B ^ qMAP are j egs p ro bable explanations of the data. The full posterior 
probability density function p{9 |D,Z) is used to generate summaries that estimate the uncertainty 
of each parameter of the model (see Section I3.2.2P . 

The observed counts above background D in the RHESS I data for bo t h flar es are large enough 



(> 30 counts in all but the very highest energy-loss bins, IWassermanl . 120031 ) that the Poisson 
distributions in Equation [T7] can be approximated by Normal distributions with mean and variance 
both equal to Ci(9 B ). Therefore, the logarithm of the posterior is approximately 



J* (Di-Ci(d B )Y 
lnp(^|D,X) oc £ l ' '^ >> (18) 



where % < n is the number of energy loss bins at which the number of counts is large enough 
that the Gaussian approximation is valid. Therefore the hypersurface formed by the Bayesian 
posterior probability density function is closely related to the x 2 -hypersurface of Equation [TUl To 
estimate 9 M and the less probable explanations of the data we turn to Markov chain Monte 
Carlo methods to efficiently explore the posterior probability density function. Note that the full 
posterior assuming the Poisson likelihood Equation 1171 was used in the analysis, and not Equation 
T8| since Equation [17] is more appropriate and the Markov chain Monte Carlo method applied 
to Bayesian data analysis does not require Normal distributions in order to generate uncertainty 
estimates. 
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We note that a similar application of Bayesian data analysis techniques was implemented 
to generate values and uncertainty estimates in the recovery of the differential emission measure 



(DEM) from emission line spectra. iKashvap fc Drakd (|1998l ) recast the DEM recovery problem 



using Bayes' theorem and modeled the full DEM as a set of emissivities and elemental abundances 
in a fixed number of temperature bins. This model is convolved with the contribution functions 
of the emission lines observed to generate a predicted emission. The parameter space describing 
the DEM is explored using a Markov chain Monte Carlo technique. The advantage of the Bayesian 
data analysis approach in DEM recovery is that it provides confidence limits on the most probable 
DEM at each temperature, thus allowing a determination of the significance of apparent structures 
that may be found in a typical reconstruction. 



3.2.1. Markov chain Monte Carlo methods for posterior sampling 

Having written down the posterior, the remaining step in the calculation is to sample from the 
posterior and calculate posterior probabilities. A brute force calculation of posterior probabilities 
can be prohibitively computationally expensive in medium or high dimensional spaces. For example, 
explicitly calculating the posterior probability density using ten different values in each of the 
seven parameters for either of the two flare models used here would require 10 evaluations of the 
posterior function. We adopt a more practical approach by using a Markov chain Monte Carlo 
method to find samples from the posterior probability density function. MCMC methods allow 
for the efficient mapping of Bayesian posterior probability density functions in multi-dimensional 
parameter space. After some initial period (known as "burn-in"), the Markov chain returns samples 
directly proportional to their probability density as defined by the Bayesian posterior, that is, 
the equili brium distribu tion of the Markov chain is the same as the posterior probability density 



function ({Gregory! 120051 ). In general, it is desirable for the Markov chain to have "rapid mixing", 
that is, it quickly reaches its equilibrium distribution. Many different MCMC algorithms have been 
designed in order to achieve rapid mixing. In this paper, we implement a parallel tempering MCMC 
algorithm (see Appendix [A] for more details). Tabled] show the priors used for each variable and 
the range of values of 9 for each flare. Assessing when the post burn-in state has been achieved 
can be found by examining the s amples. In this paper, the Gelman R diagnostic is used to assess 



convergence (JGelman et al.l 120031 . see Appendix 



3.2.2. Summaries of the posterior probability density function 

The Bayesian/MCMC summary probability density functions for a single parameter 9i in the 
set 9, (1 < i < Ng) are found by integrating the posterior probability distribution (Eq. [TBI) over all 
the other variables, i.e., 

p(9i)= [ P (H\T>,l)d9 1 ...d9 i - 1 d9 i+1 ...d9 Ne . (19) 
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This distribution is called a marginal distribution, and it is the probability density function for the 
variable Oi given all the likely values of all the other variables. The marginal distribution is used to 
calculate uncertainty estimates to 6{. Values to the 68% and 95% uncertainty are calculated using 
the definition of the uncertainty interval given by Equation 1141 with the function F{0i) substituted 
with the marginal distribution p(0i). The uncertainties quoted for this method in Tables [2] and [3] 
are given as 

9% - median [p(0;)] ,6 H \ q - median [p(9i)] (20) 

where 9 L, 9 H \ q are defined using Equation Q3] (substituting p(9i) for F(9i)) for q = 0.68 and 
q = 0.95 and median [p(9i)] is the median value of the marginal probability density function p(6i). 
Note that this definition of the interval does not necessarily include the mean or the mode. 



4. Results 

4.1. 19 January 2005 

Figures [3l [U [5] and Table [2] show the results for each of the four uncertainty estimation methods 
under consideration using the data and electron spectral model for the 19 January 2005 flare, as 
described in Section [2j Figures EJ H] and Table [2] show that the difference between the 9 and 9 

values are much less than the 68% uncertainty estimates. For each variable, the lower and upper 68% 
(and 95%) uncertainty estimates found by each uncertainty estimation method have approximately 
the same magnitude. Comparing across methods, it can be seen that each also gives approximately 
the same uncertainly estimates. The ratio of the 95% uncertainty estimate to the 68% uncertainty 
estimate are all close to 1.96, as expected from distributions of measurements which are close to 
being Normally distributed. In addition, Q-Q plots of all seven marginal distributions obtained 
from the Bayesian analysis (see Appendix [C]) show that each of them is approximately Normally 
distributed. 

Figure [5] plots two-dimensional marginal distributions arising from the Bayesian/MCMC anal- 
ysis for every pair of parameters in the spectral model (the priors used in the Bayesian/MCMC 
approach can be found in Table |4|). It shows the effect each parameter has on the value of the other 
when finding highly probable parameter values to 9. Next to each plot the Spearman rank corre- 
lation coefficient for the indicated variables is shown. It can be seen that all the two-dimensional 
marginal distributions are elliptical, and the majority of them show that the probability of getting 
a particular parameter value is weakly correlated with the value of any other parameter. The 
exceptions to this for this flare are the emission measure (EM) and plasma temperature (kT) de- 
pendency, the dependency of the spectral normalization Fq on the low-energy cutoff E c and the 
power law index S±, and the E c versus 5± correlation. 

The first of these dependencies is anticipated through the definition of the thermal emission of 
the plasma (Equation [7]) , and the second two arise from the definition of the normalization. The 
normalization factor Fq for this flare is defined as the total integrated electron flux over all energies, 
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and therefore clearly depends on the values of E c and 5\ (see Section [2]). Figure [5] also shows a 
correlation between E c and Si. This is obtained because the rate at which the X-ray spectrum 
flattens below E c depends on the value of Si. The spectrum flattens more rapidly with decreasing 
photon energy for a steeper electron distribution (larger Si) than for a flatter electron distribution. 
Therefore, for a given X-ray spectrum, a larger value of Si requires a higher value of E c to obtain 
the best fit to the spectrum. A similar correlation, for the same reason, is found between E^ and 
82 in the fit to the July 23 flare spectrum (Figured]). 

Figure[6|a) shows the (scaled) electron flux spectrum as a function of energy for the Bayesian/MCMC 
analysis. Figure Efb) shows the ratio of the best fit electron spectrum to the 68% and 95% uncer- 
tainty estimates. Figures UK a) and (He) show the probability density functions for total integrated 
electron number flux and electron power derived from the Bayesian/MCMC results. Uncertainty 
estimates for the electron flux spectrum as a function of energy are found in the following way. 
The electron flux spectrum for each Bayesian/MCMC-derived sample is calculated. The spectra 
are then ranked according to their posterior probability. The 68% curves are found by finding the 
highest and lowest values to the electron flux spectrum in each energy bin for the top 68% most 
probable samples (the 95% curves are found similarly), yielding the uncertainty estimates as shown 
in Figure ED(a). In each energy bin, the upper and lower uncertainties are approximately symmetric 
around the best (9 M ) value. Further, the probability density functions for the electron number 
flux and power (Figures IHa) and (He)) are also approximately symmetrical around the mean and 
mode. This is not too surprising since the probability density functions (Figures [3l H|) for each 
parameter in the fit are also approximately symmetrical. Finally, the uncertainties in the values to 
the electron number and power are also well constrained. 



4.2. 23 July 2002 

Figures El [HJ [9] and Table [3] show the results for each of the four uncertainty estimation meth- 
ods under consideration using the data and electron spectral model for the 23 July 2002 flare, as 
described in Section [2j It is clear from Figures [8] and [9] that the x 2 -hypersurface (or equiva- 
lently, the Bayesian posterior hypersurface - see Section 13.2. ip with respect to this model is quite 
different from that seen in the 19 January 2005 flare (Figures 2] and [5|) . The mode values in the 
Bayesian/MCMC marginal distributions are noticeably shifted with respect to the Monte Carlo 
distributions. This is because the Bayesian/MCMC marginal distributions in Figures [7] and [8] are 
formed by integrating over a structured seven-dimensional space (Figure [9]). The mode of the 
one-dimensional marginal distributions need not be at the Q MAP or 9 value. Note however from 
Table [3] that the MAP value is close to the 6 value, which is to be expected given the priors used 
in setting up the Bayesian posterior (see Appendix [B]) and the close correspondence between the 
X 2 -hypersurface (Equation I10p and the Bayesian posterior (Equation I18p . 

Figures [7] (thermal model parameters) and [8] (non-thermal model parameters) show that the 
uncertainty estimates for specific parameters can depend on the uncertainty estimation method 
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Fig. 3. — Results from each of the four uncertainty analysis methods (Section [3]) for the parameters 
of the thermal component of the total emission (a) EM and (b) kT, from the model fit to the 19 
January 2005 flare data. The dashed curve is the value of x 2 found by the % 2 -mapping method 
(values are indicated by the right-hand plot axis) . The normalized frequency distribution of values 
found by the Monte Carlo method is shown as a histogram (dot-dashed line). The marginal 
probability density function arising from the Bayesian/MCMC method is shown as a histogram 
(solid line). Values to these histograms are indicated by the left-hand plot axis. The horizontal 
lines show the uncertainty estimates calculated via the methods indicated (from top to bottom 
- covariance matrix, Bayesian/MCMC, Monte Carlo, and x 2_ma PPi n g)) with the 



and 95% 

uncertainty estimates indicated by larger and smaller vertical lines that cross those lines. The 
best-fit value 6 found via nonlinear least-squares minimization (Section 13. ip is indicated by square 
plot symbols. The MAP value Q MAP is indicated by a x -symbol. The mean, mode and median 
values calculated for each of the two distributions (arising from the Bayesian/MCMC and Monte 
Carlo analyses) are indicated by asterisks, diamonds and triangles respectively. These symbols are 
separated vertically scattered for clarity. 
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(a) F , 2005 January 19 



b) 6 , 2005 January 19 
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Fig. 4. — Results from each of the four uncertainty analysis methods (Section [3]) for the parameters 
of the nonthermal component of the total flare emission (a) Fq, (b) 8± and (c) E c (see Eq. [9]) from 
the model fit to 19 January 2005 flare data. The type of data plotted, plot symbols and lines have 
the same meaning as in Figure El 
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Fig. 5. — Two dimensional marginal probability density functions for the parameters of the model 
used to fit the spectrum of the 19 January 2005 flare. These plots are found by integrating the 
posterior probability density function (found by the Bayesian/MCMC algorithm) over all the pa- 
rameters excepting those indicated on the x- and y-axes. This is the extension into two dimensions 
of the definition of the one-dimensional marginal distribution function given by Equation [19] in 
Section 13.2.21 Each of these plots in this figure shows how the posterior probability density of the 
value of a given parameter depends on the value of another parameter, and so help visualize the 
shape of the full posterior probability density function. Indicated parameter ranges are the low- 
est and highest values found by the Bayesian/MCMC algorithm. Darker tones indicate a greater 
probability density. The number on the upper right of each plot is the Spearman rank correlation 
coefficient for the two parameters. For the 19 January 2005 flare, the distributions are all approxi- 
mately elliptical. The majority of the distributions are weakly correlated; a minority (EM versus 
kT, and Fq, 5\ versus E c , Fq versus 5\) show a high degree of correlation. The reasons for these 
strong correlations are discussed in Section [4.11 
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Table 2. Parameter values and uncertainties derived for the four uncertainty estimation methods 

applied to the 19 January 2005 flare spectrum, as described in Section [2j The final column 

"Ratio" is defined as the ratio of the ±95% uncertainties to the ±68% uncertainties; for an exact 

Normal distribution the entry in this column would be 1.96, 1.96. Two ratios are quoted in order 

to reveal the presence of any relative asymmetry in the upper and lower uncertainty estimates, if 

present. See Section [3] for a detailed description of how the uncertainty estimates are found for 

each method. 



Parameter 


Method 


Value a 


Uncertainties 


Ratio 








68% 


95% 




EM (10 49 cm- 3 ) 


covariancc matrix b 


2.31 


±0.14 


not calculated 


not calculated 




X 2 -mapping c 


" 


-0.14, +0.15 


-0.27, +0.31 


1.94, 2.05 




Monte Carlo d 


" 


-0.17, +0.12 


-0.30, +0.27 


1.75, 2.31 




Bayesian/MCMC° 


2.30 


-0.14, +0.15 


-0.27, +0.31 


1.96, 2.04 


kT (kcV) 


covariance matrix 


2.03 


±0.02 


not calculated 


not calculated 




X 2 -mapping 


" 


±0.02 


±0.04 


1.99, 2.01 




Monte Carlo 


" 


±0.02 


-0.03, +0.04 


2.13, 1.84 




Bayesian/MCMC 


2.03 


±0.02 


±0.04 


1.98, 2.03 


F 


covariance matrix 


0.17 


±0.01 


not calculated 


not calculated 


(total integrated electron flux 


X 2 -mapping 


" 


±0.01 


±0.02 


1.90, 2.10 


10 35 electrons sec - 1 ) 


Monte Carlo 


" 


±0.01 


±0.01 


1.87, 2.17 




Bayesian/MCMC 


0.16 


±0.01 


±0.02 


1.94, 1.96 


h 


covariance matrix 


3.57 


±0.03 


not calculated 


not calculated 




X 2 -mapping 


" 


±0.04 


-0.07, +0.08 


1.95, 2.04 




Monte Carlo 


" 


-0.02, +0.04 


-0.05, +0.07 


2.15, 1.89 




Bayesian/MCMC 


3.58 


-0.03, + 0.04 


-0.06, +0.07 


1.88, 2.04 


E c (keV) 


covariancc matrix 


105 


±3 


not calculated 


not calculated 




X 2 -mapping 


" 


±1 


±8 


2.00, 2.00 




Monte Carlo 


" 


-3, 4 


-6 , +7 


2.10, 1.91 




Bayesian/MCMC 


107 


±4 


-7, +8 


1.85, 2.00 



a The covariance matrix, Monte Carlo and x 2 -nrapping methods all start from the same parameter value 9 where 
X 2 is minimized. For the Bayesian/MCMC approach, the "maximum a posteriori" value Q MAP is quoted. 

b See Section 13.1.11 and Equation 1111 for the definition of the parameter uncertainty for the covariance matrix 
method. 

c See Section 13.1.31 and Equation [15] for the definition of the parameter uncertainty for the x 2 -mapping method. 

d See Scction l3. 1.21 and Equation [13] for the definition of the parameter uncertainty for the Monte Carlo method. 



c See Section 13.2.11 and Equation 1201 for the definition of the parameter uncertainty for the Bayesian/MCMC 
method. 
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Fig. 6. — Electron spectrum results for the flare- injected electrons arising from the 
Bayesian/MCMC method for the 19 January 2005 flare, (a) Electron spectrum (flux (in units 
of erg keV~ s _1 ) multiplied by I? 3,58 , where 3.58 is the MAP estimate <5i, the power law index 
of the flare-injected electron spectrum - see Table [2]) with 68% and 95% credible interval spectra 
indicated by the dashed and dotted lines, respectively. The electron flux spectrum corresponding to 
qMAP j g i nc ii ca t ec i by the solid line, (b) 68% and 95% credible intervals relative to the # elec- 

tron flux spectrum. In plots (a) and (b) curves with negative gradients indicate a behavior steeper 
than E~ 1 and positive gradients indicate a behavior shallower than E~ 1 . Note also that the MAP 
spectrum extends to its low energy cutoff value; other lower probability spectra extend to values of 
E c , which may be different to the MAP value of E c . (c) Flare injected electron power probability 
density function, with 68% and 95% credible intervals indicated; the distribution mean/mode is 
indicated by the solid/dot-dashed vertical line. The total integrated electron flux injected by the 
flare is given in Figure IHc). 
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used. The methods used are influencing the uncertainty estimates for some parameters (Table [3]). 
These uncertainty estimates behave quite differently from those expected from a Normal distri- 
bution, with the ratios of the 95% to 68% uncertainty estimates very different from 1.96. The 
reason for this is apparent when considering the two-dimensional Bayesian posterior marginal dis- 
tributions as shown in Figure Many of the distributions are structured, asymmetric, and show 
extended tails compared to those derived from the hypersurface of the 19 January 2005 analysis. 
The low-energy cutoff in particular shows significant deviation from a simple Normal distribution, 
as does the break energy E\, and the slope of the spectrum above the break energy, parameterized 
by 82- Many pairs of parameters have high magnitude correlation coefficients indicating strong 
interdependence of one value on another. Further, note that the correlation of E c with all other 
parameters is relatively weak. This indicates the relative independence of the low-energy cutoff 
from other features in the model, given the data. 

Figure EJe) and [9] show that below around 25 keV, all values of E c are approximately equally 
likely, but also that E c < 25 keV does not constrain likely values of the emission measure EM, the 
thermal temperature kT, the normalization A and the lower power-law index 6\. This leads to a 
wide range of possible electron-flux spectra at lower energies, the effect of which leads to wide 68% 
and 95% credible intervals of Figure [TUlf a). The uncertainty estimates for the electron flux in Figure 
[6]^ a) also show a widening at lower energies, but it is much less pronounced compared to that in 
Figure [TUlf a) . The reason for this is that at lower values of E c , the other parameter values in the 
model are constrained, and so there is a restricted range of electron flux spectra that is generated. 

Figure [101 shows the (scaled) electron flux energy spectrum as a function of energy, along with 
probability density functions for total integrated electron number flux and electron power derived 
from the Bayesian/MCMC results. The wide 68% and 95% credible intervals of Figures [TOlfa. 
b) show that the electron spectrum becomes poorly constrained at low energies. Figures [TUtc) 
and (d) are the electron number and power probability density functions, respectively (found by 
integrating the flare spectrum electron flux spectrum from E c to E^). Both are asymmetric and 
show more pronounced tails when compared to the corresponding plots for the 19 January 2005 
data (Figures Hl^a) and[6^c)). This is due to the asymmetric low-energy cutoff probability density 
function which leads to a tail extending to high values in the probability density function of the 
electron number flux. Uncertainty estimates for the total number of flare-accelerated electrons and 
their energy are given in Figures fTOT c. d). The probability density function for the energy can 
be integrated to determine lower limits to the energy contained in the flare-accelerated electrons 
whilst simultaneously supplying a probability estimate. The cumulative probability distribution 
function for the energy shows that there is a 95% probability that the energy in the flare-accelerated 
electrons is greater than 10 28,0 erg sec -1 , and a 68% probability that it is greater than 10 28 ' 2 . 



27 



Table 3. Parameter values and uncertainty estimates derived for the four uncertainty estimation 
methods applied to the 23 July 2002 flare spectrum, as described in Section [2l The final column 
"Ratio" is defined as the ratio of the ±95% uncertainties to the ±68% uncertainties; for an exact 
Normal distribution the entry in this column would be 1.96, 1.96. Two ratios are quoted in order 
to reveal the presence of any relative asymmetry in the upper and lower uncertainty estimates, if 
present. See Section [3] for a detailed description of how the uncertainty estimates are found for 
each method. The entry 'not determined" indicates that the value was not determinable by the 

method. 



Parameter 


Method 


Value a 


Uncertainties 


Ratio 








68% 


95% 




EM (10 49 cm- 3 ) 


covariance matrix b 


2.16 


±0.08 


not calculated 


not calculated 




X 2 -mapping c 




±0.04 


±0.08 


2.05, 1.99 




Monte Carlo d 




-0.05, 0.03 


-0.09, 0.07 


1.82, 2.28 




Bayesian/MCMC e 


2.17 


±0.04 


±0.08 


1.89, 1.96 


kT (keV) 


covariance matrix 


3.18 


±0.03 


not calculated 


not calculated 




X 2 -mapping 




±0.01 


±0.02 


1.97, 2.13 




Monte Carlo 




±0.01 


-0.02, 0.03 


2.20, 1.87 




Bayesian/MCMC 


3.18 


±0.01 


±0.03 


1.93, 1.92 


A 


covariance matrix 


0.028 


±0.004 


not calculated 


not calculated 


(electron flux at 50 keV, 


X 2 -mapping 




-0.003, 0.002 


-0.006, 0.005 


2.15, 1.94 


10 35 electrons (sec keV)" 1 ) 


Monte Carlo 




-0.002, 0.003 


-0.005, 0.005 


2.21, 1.74 




Bayesian/MCMC 


0.028 


-0.003, 0.002 


-0.006, 0.004 


2.09, 1.76 


Si 


covariance matrix 


3.40 


±0.16 


not calculated 


not calculated 




X 2 -mapping 




-0.14, 0.10 


-0.36, 0.17 


2.61, 1.78 




Monte Carlo 




-0.14, 0.12 


-0.34, 0.19 


2.52, 1.61 




Bayesian/MCMC 


3.41 


-0.13, 0.08 


-0.33, 0.13 


2.55, 1.55 


E b (keV) 


covariance matrix 


256 


±135 


not calculated 


not calculated 




X 2 -mapping 




-77, 147 


-123, 686 


1.59, 6.67 




Monte Carlo 




-77, 253 


-121, 1319 


1.58, 5.22 




Bayesian/MCMC 


269 


-147, 5615 


-217, 1239 


1.47, 2.01 


,h 


covariance matrix 


3.92 


±0.11 


not calculated 


not calculated 




X 2 -mapping 




-0.08, 0.13 


-0.13, 0.78 


1.67, 5.67 




Monte Carlo 




-0.07,0.23 


-0.12, 3.27 


1.74, 14.2 




Bayesian/MCMC 


3.93 


-0.11, 0.58 


-0.18, 1.92 


1.57, 3.33 


E c (keV) 


covariance matrix 


32.0 


±24.091 


not calculated 


not calculated 




X 2 -mapping 




-5.78, 5.05 


not determined, 12.1 


not determined, 2.4 




Monte Carlo 




-6.86, 7.37 


-20.7, 15.9 


3.02, 2.16 




Bayesian/MCMC 


31.2 


-16.1, 11.7 


-23.1, 19.1 


1.44, 1.63 



a The 'covariance matrix', 'Monte Carlo' and 'x -mapping' methods all start from the same parameter value 8 where \ 
is minimized. For the Bayesian/MCMC approach, the "maximum a posteriori" 8 M value is quoted. 
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b See Section |3. 1 , 21 and Equation 1131 for the definition of the parameter uncertainty for the covariance matrix method. 
c See Section 13,1.11 and Equation [TT] for the definition of the parameter uncertainty for the x 2 -mapping method. 
d See Section |3. 1 .31 and Equation 1151 for the definition of the parameter uncertainty for the Monte Carlo method. 
e See Section |3. 2. II and Equation 1201 for the definition of the parameter uncertainty for the Baycsian/MCMC method. 
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As was noted in Section 12.21 a different spectral normalization was used in the analysis of the 
23 July 2002 flare compared to the 19 January 2005 flare. The package OSPEX implements the 
spectral normalization of the 19 January 2005 model spectrum using the integrated normalization 
factor, Fq = AEq~ f[8\ — 1)- This implementation of the flare spectral model therefore introduces 
a parameter dependence into the % 2 -hypersurface between the normalization A, the low-energy 
cutoff and the spectral index S±. However, since the low-energy cutoff for the 19 January 2005 flare 
is relatively well defined, the integrated flux Fq is relatively well defined, and the MCMC algorithm 
can explore the x 2 -hypersurface as a function of Fq and E c with no difficulty. However, the low- 
energy cutoff is not well defined for the 23 July 2002 flare, and so the range of values to Fq is large. 
Therefore when using the implementation of Equation [9] used in the analysis of the 19 January 2005 
flare, the parameter space that must be covered by the MCMC algorithm is large due to the inherent 
dependence of Fq on E c . This was found to be prohibitive to an efficient MCMC search, and so an 
alternate implementation of Equation [9] was created for OSPEX (re-paramet erization of the fitti ng 



function is a recommended tactic in creating better search spaces for MCMC (jGelman et al.ll2003l )). 
In this implementation, the normalization factor used to describe the spectrum is A, the value of 
the spectrum at the pivot value E p . Moving to a different hypersurface for the same problem 
greatly improved the efficiency of the MCMC algorithm. 



5. Discussion 

5.1. Comparison of uncertainty analyses 

The uncertainty analyses performed on both data-sets shows that the shape of the x 2 -hypersurface 
has a significant effect on the values of the uncertainties found. All the uncertainty estimates found 
for the spectral parameters describing the 19 January 2005 flare data are similar, regardless of the 
method. The uncertainty estimates found for the spectral parameters describing the 23 July 2002 
flare data depend on the method chosen. 

Since the data have a large number of counts at almost all energies, the hypersurfaces described 
by Equation [10] and Equation [18] are almost identical. The two-dimensional marginal distributions 
for the 23 July 2002 flare data (Figure [9]) shows structures which are not simple two-dimensional 
Normal distributions, and, since the two hypersurfaces described by Equation [18] and Equation [TUl 
are almost identical, the x 2 -hypersurface must have structures which are not simple two-dimensional 
Normal distributions. This means that one or more of the assumptions that lead to the assertion 
that the probability distribution for 50 o b s is a multivariate Normal distribution around 9 does not 
hold for this model applied to these flare data (Section l3.1.ip . The non- Normal distribution shapes 
of Figure [9] suggest that the assumption that the spectral model is li near (or at least l ocally so 



within the range of the desired uncertainty calculation) is not satisfied ([Press et al.l Il992l . p. 690). 
Hence, the covariance matrix and x 2 " ma PPi n g methods cannot be expected to give reliable and 
consistent estimates in this case. 
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The shape of the x 2 -hypersurface also influences the results of the Monte Carlo method. This 
can be seen in the results for the low-energy cutoff in the 23 July 2002 data-set (Figure [8b). It is 
expected that below a given energy E p i ateau , all values of the low-energy cutoff are equally likely. 
This is because in this energy range the number of counts due to thermal emission greatly exceed 
the number due to the flare- injected electron flux spectrum, and so changing one value of E c over 
another makes no difference to the fit to the data - the value of x 2 > or equivalently, the Bayesian 
posterior probability, are unaffected. Therefore, all values below E p i a t ea u are equally likeljUI. The 
Monte Carlo method results do not show this; the results are clustered around the best-fit value and 
do not show the extension to lower energies as expected. Hence the uncertainty estimate arising 
from the Monte Carlo method does not conform to our prior expectation of what it should report. 

In contrast, the Bayesian posterior hypersurface for the 19 January 2005 shows simple Normal- 
like one-dimensional distributions (and so the assumptions behind the covariance matrix and \ 2 - 
mapping methods are approximately true) and give similar answers. The Monte Carlo method 
(Section 13.1. 3D relies on finding local minima to simulated data which is statistically similar to 
the original data. This method works well in the 19 January 2005 analysis as the shape of the 
hypersurface (Figure [5j) is dominated by a nearly Normal single minimum, a feature the method 
repeatedly finds in all the similar x 2 -hypersurface. The x 2 - ma PPi n g method does agree with the 
Bayesian/MCMC result in that the x 2 - ma pping method does indicate that below a certain value 
(Epi a t e au), all values of the low-energy cutoff are equally likely. However, the method cannot give a 
lower limit to the 95% uncertainty estimate since at no point does 5\ 2 = 4 for E < E™ in (Section 
15X21) . 

The Bayesian/MCMC method samples the parameter space via the posterior probability and 
the Markov chain Monte Carlo algorithm (Section 13.2. lj) . The Bayesian interpretation of the 
posterior probability means that the parameter samples are found in proportion to how well they 
describe the data (values of 9 that have lower probability are less likely explanations of the data). 
The method does not make any assumptions about the nature of the hypersurface, as the other 
three methods do. Hence it agrees with the results from the methods of Section 13.11 when applied 
to simple hypersurfaces where the assumptions made by those methods are valid, but generates 
different results when those assumptions do not hold. Therefore, the Bayesian/MCMC method 
can, in principle, be used without having to invoke any special knowledge of the shape of the 
hypersurface and without making some simplifying assumptions. 
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5.2. Probability density functions of the parameters of the 23 July 2002 electron 

spectrum model 

Figures [7] and [8] show the marginal probability density functions of the parameter values arising 
from a Bayesian/MCMC treatment of the data analysis problem. It is notable that the distribu- 
tions for Eb, 62 and E c are distinctly different from more symmetrical and Normal distribution-like 
distributions of the other parameters in the fit. The break energy Ej, and the power law index 
above the break 62 are highly correlated (Figure [9]) over a wide range of values. As Ef, increases, 
the value of 82 increases. The mild curvature of the spectrum implied by these probability density 
functions is consistent with a wide range of near power-law electron flux spectrum models, leading 
to an ill-defined value for Ef, and softer power-law indices at higher values of Ef,. A count spectrum 
that appears to come from emission that is mildly curved with respect to the radiation from the 
thick-target interaction of a flare-injected electron flu x spectrum with a power law distribution 



could arise from an inaccurate X-ray albedo correction (IKontar et al.ll2006r ) or from a non-uniform 



ionization within the target plasma ()Su et al.l 12009 ; IKontar et al.ll2002l ) 



The low-energy cutoff also has an interesting probability density function (also reproduced by 
the x 2_ma PPi n g analysis, Figure [5fe). There is a peak in the Bayesian/MCMC low-energy cutoff 
probability density function at 31 keV, and a tail at lower energies where the thermal emission 
of the plasma dominates over the emission due to the flare- injected electron flux. We wish to 
estimate how much more likely the low-energy cutoff is close to the peak, compared to other parts 
of the probability density function. An estimate can be generated using the following procedure. If 
the probability density function of the low-energy cutoff were a Normal distribution N(E cu toff,&) 
(where N(a,b) is a Normal distribution centered at a with standard deviation b), then the total 
probability that E c lies in the range E cu toff — &, E cu toff+v is about 68%. The maximum probability 
that E c lies in a 2a wide range of values that does not overlap with the range E cu t ff — a, E cu t ff + 
is about 16%. Therefore the value of E c is about 4 times more likely to be in the range E cu toff — 
a, E cu toff + o than in a la wide range of values that does not overlap with the range E cuto ff — 
a, Ecutoff + G - Fitting the peak of the probability density function of Figure [8^e) with a Normal 
distribution yields a width a of about 5 keV. Applying the estimation procedure above on the 
probability density function of Figure [8]^e) with a = 5 keV, it is found that E c is about 1.3 times 
more likely to be in the range 25-35 keV than in any other continuous window of values 10 keV 
wide. This is weak evidence for a peak in the range 25-35 keV. 

Therefore, the probability density function is interpreted as providing evidence for the existence 
of an observable low-energy cutoff just above the region where the thermal emission dominates. If 
the low-energy cutoff was at higher energies, then the probability density function for E c would 
resemble more closely the probability density function seen in Figure [U^c) for the January 19 flare 
and therefore lower possible values to E c would lead to lower posterior probabilities (worse fits). 



7 Epi a t e au can also be interpreted as the energy below which no further information is available that can be used 
to better constrain a lower limit to the low-energy cutoff. 
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If the low-energy cutoff was present at energies where the thermal emission dominates, then no 
peak in the probability density function for E c would be seen. Lower values would account for 
more of the flare-injected spectrum, and so lower values would be more probable. The probability 
p{E c ) would eventually plateau at some energy E v i ateau since the emission due to the flare- injected 
electron flux would be far less than the emission due to the thermal plasma below E p i ateau , making 
all values of E c equally likely, as there is nothing to distinguish one value from another. However 
the observed p(E c ) is a combination of both; a peak in the probability density function with an 
approximately constant probability density at lower energies. 



5.3. Flare electron number and energy probability density functions 

The Bayesian/MCMC method allows for the construction of probability density functions for 
each flare (Figures Ufa) and[6jc), [TOlfc. d) of the number of flare-accelerated electrons and the 
energy they carry, fully expressing the correlated dependence of one variable on another (Figures 
[5l [9]). Since the result is another probability density function, credible intervals for the number of 
electrons and their energy can also be calculated. In contrast, taking the set of 68% upper model 
parameter uncertainty estimates (or the other model parameter uncertainty estimates) from the the 
other methods cannot be used to calculate the corresponding 68% upper uncertainty estimate for 
the number of electrons and their energy. This is because there is no guarantee that that point on 
the x 2 -hypersurface has a significant non-zero probability (or equivalently, lies in a hightly probable 
region of the model parameter hypersurface). In relatively simple hypersurfaces this may be true, 
but in highly correlated hypersurfaces such as in the analysis of the 23 July 2002 flare presented 
here, it may not be. As far as we are aware, this is the first time that flare electron number and 
energy probability density functions have been estimated from data. 

A significant difference between the two flares studied is the uncertainty with which the model 
parameters are known. This leads to significant differences in how well the gross properties of the 
flare are known. The low-energy cutoff is not well constrained for the 23 July 2002 flare, leading 
to 68% and 95% credible intervals in the flare electron number and energy probability density 
functions that span orders of magnitude. Notably, the 23 July 2002 probability density functions 
are highly asymmetric and so lower values of flare electron number and energies are much less likely 
than higher values. It is interesting to note that there is a peak in the energy probability density 
function for the 23 July 2002 flare, even although there is a non-zero probability for E c down to 
the lower limit given by the prior for the low-energy cutoff. This is due to the peak in the marginal 
probability density function of E c , which therefore defines a more probable total flare energy than 
those arising from the lower probability range E c < E p i ateau . 

The estimate of the actual number of electrons and the energy they carry is also dependent 
on systematic errors related to the calibration of each of RHESSI detectors with each other. As 
was noted above, the systematic errors in the individu al PHA bins are small compared to th e 



systematic error in the overall sensitivity of each detector (JMilligan &; Dennidl2009l : ISu et al. 



2011 
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This means that the shape of the flare-accelerated electron spectrum suffers from a smaller error 
compared to the integral under the curve of the flare-accelerated spectrum. We therefore expect 
that the broad qualities of the shapes of the flare electron number and energy distributions will 
remain unchanged for each of the two flared studied; the 19 January 2005 results will remain 
approximately symmetric, and the 23 July 2002 results will remain quite asymmetric. We estimate 
that allowing for a 10% - 30% error in knowledge of the sensitivity of each detector would smooth 
out the distribution peak, and add another 0.1-0.2 in the logarithm (approximately) of the widths 
of the probability density functions. This estimated uncertainty is substantially more than the 95% 
estimated uncertainty in the case of the 19 January 2005 flare, but is substantially less than the 
95% estimated uncertainty for the 23 July 2002 flare. This suggests that the uncertainty in the true 
value of the low-energy cutoff is a more important limiting factor in understanding the electron 
and energy content in RHESSI-observed flares than the detector calibration uncertainty. 



5.4. Expanding the analysis 

It is common in RHESSI data analysis to remove a background component from the observed 
count data to yield an estimate of the counts due solely to the flare. This background-subtracted 
data is then used in further analysis. Strictly, models for the background and the flare should be fit 
simultaneously since the observed counts are due to the background and the flare simultaneously. 
Therefore, the first improvement we will make is to fit both the flare response and background 
simultaneously. This will be done by including a simple parameterization of the pre- and post-flare 
hard X-ray flux observed by RHESSI into the flare model. The parameters of the background 
model will also require their own priors. The inclusion of a background model in the fit is expected 
to have an effect an higher energies, where the signal-to-noise ratio of the flare-accelerated electrons 
are smaller, such as in smaller flares. 

The analyses presented here made use of data from one single detector. Our second improve- 
ment to the existing analysis will be to including data from more than one detector, which will 
increase the signal to noise ratio. In order to use data from more than one detector, information 
about the relative calibration of each detector will have to be included. This will be incorporated 
into priors for each detector that express the degree of uncertainty in their calibration. Since each 
detector is observing the same flare, the flare model will be the same across detectors. The poste- 
rior will be a product of the priors for the flare model plus background, a likelihood function for 
each detector, and a prior function expressing the degree of uncertainty in their calibration. The 
resulting posterior will express the increased knowledge that comes with a larger number of counts, 
but also the uncertainty in their relative calibration. 

We note also that Bayesian data analysis provides a framework that can be used to compare 
the explanatory power of differ ent models of the data whilst taking into account the number and 



type of variables in each model (|GregoryM2005l ). We will use Bayesian model comparison techniques 



to determine if RHESSI data can distinguish between different effects that may contribute to the 
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observed spectra. In particular, we will re-analyze the 23 July 2002 data prese nted here 



a model that incor porates the non-uniform ionization of the thick-target plasma (|Su et al.l 12009 ; 



using 



Kontar et al.l 120031 ). Such a model produces a curvature in flare-accelerated electron spectrum 
which may explain the high correlation between the break energy Eb and the value of 82 (Section 



6. Conclusions 

This paper describes in some detail four methods that can be used to estimate the uncertainties 
in parameters of flare models fit to RHESSI hard X-ray flare data. Three of the four methods - 
covariance matrix, Monte Carlo, and x 2 -mapping - measure scale-sizes in the x 2_ hypersurface (or 
related hyper surf aces) and call them uncertainty estimates. We have shown that care must be taken 
in relying upon these uncertainty measurements, as we have seen that they need not agree with our 
expectation of what an uncertainty estimate should report, or with each other. The fourth method, 
Bayesian data analysis, can answer the question "what is the uncertainty in this parameter?" by 
calculating a probability density function for that parameter through the marginalization procedure 
of Section 13.2.21 without making any further assumptions about the number of counts in each bin 
(see Section l3.2.ip . The fourth method broadly agrees with the other three in the case of the 19 
January 2005 flare. Each method generates different uncertainty estimates for the 23 July 2002 
flare. 

The source of the different uncertainty estimates is the shape of the x 2 -hypersurface parame- 
terized by the flare model, hypersurfaces that broadly conform to the assumptions underlying the 
covariance matrix, Monte Carlo, and x 2 -mapping methods yield consistent uncertainty estimates 
that agree with each other and those from the Bayesian/MCMC approach. Conversely, hyper- 
surfaces that break those assumptions yield method-dependent results. The Bayesian/MCMC 
approach makes no assumptions on the nature of the hypersurface. Further, the position of the 
low-energy cutoff in relation to the region where thermal X-ray emission dominates is crucial in 
determining the shape of the hypersurface. Most flares are thought to have a low-energy cutoff 
close to or at the region of thermal emission dominance. The Bayesian/MCMC method presented 
here handles both flare analyses without regard to the location of the low-energy cutoff, and makes 
no assumption about the x 2_ hypersurface or Bayesian posterior probability hypersurface. The 
Bayesian/MCMC method was the only method to generate an uncertainty estimate of the low- 
energy cutoff that reflects our intuition of how it is constrained by the data, for both flares studied. 
Since the x 2 - m apping approach does partially map the space around 9, it is perhaps the best of the 
three non-Bayesian based methods that can give an indication that the x 2_ hypersurface contains 
features that are not similar to Normal distribution shapes. If the x 2 -hypersurface does contain 
features not anticipated by the covariance matrix, Monte Carlo, and x 2_ma PP m g methods, then 
we suggest a Bayesian/MCMC approach is warranted if reliable uncertainty estimates are desired. 

The 23 July 2002 flare shows evidence for the existence of a low-energy cutoff in the range 
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25-35 keV, just above the region where the thermal emission dominates. The probability density 
function of the low-energy cutoff shows significant non-zero probability below 25 keV, and zero 
probability above 50 keV. This peak is important, as it leads to highly asymmetric probability 
density functions for the total number of flare electrons accelerated by the flare, and the energy 
they carry, in which the upper limit to these quantities are poorly constrained. In each of these 
quantities, the 95% upper credible limit is orders of magnitude larger than the MAP value, whilst 
the 95% lower limit is within one order of magnitude of the MAP value. In comparison, the MAP 
values for the same quantities of the 19 January 2005 flare lie are approximately centered within 
a tenth of a decade. This points to the importance of the low-energy cutoff probability density 
function in determining the quality of our knowledge of the gross properties of the flare. 

Further work will involve improving the modeling of RHESSI observations by including data 
from other RHESSI detectors, incorporating the simultaneous fitting of the background emission at 
the same time as the flare model, and testing different models of flare emission for the same flare. 

This work was supported by a NASA ROSES award made under the opportunity NNH09ZDA001N- 
SHP entitled "Investigation of the low energy cutoff in solar flares", and by the HESPE (High 
Energy Solar Physics Data in Europe) collaboration. We are grateful to D. van Dyk and C. A. 
Young for their helpful suggestions. CHIANTI is an Atomic Database Package for Spectroscopic 
Diagnostics of Astrophysical Plasmas. It is a collaborative project involving the Naval Research 
Laboratory (USA), the University of Florence (Italy), the University of Cambridge and the Ruther- 
ford Appleton Laboratory (UK). 
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A. Parallel tempering Markov chain Monte Carlo algorithm 

A significant problem in MCMC is ensuring that the posterior is explored sufficiently. The 
first MCMC algorithms used in this study did not generate the expected marginal probability 
distribution of the low-energy cutoff E c for the flare of 23 July 2002. The distribution arising from 
these MCMC algorithms showed a single peak with p(E c ) = below some value. The expected 
distribution contains a plateau region of approximately constant non-zero probability density for 
values E c < E v i ateau for some value of E v i atea u determined from the data (see also Section 15. ip . 
The difference between the expected distribution and those derived from the MCMC algorithm 
may be due to either insufficient exploration of the posterior by the MCMC algorithm, or to some 
previously unexpected feature in the flare spectrum. To test these explanations, a new MCMC 
algorithm was implemented to more fully explore the parameter space of the posterior distribution. 

The parallel tempering algorithm allows one to explore the parameter space by optionally 



making easier moves in related spaces (JGregoryl 120051 ) . Parallel tempering is based on simulated 



tempering. This scheme mimics the physical process of annealing, whereby a metal is heated 
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and cooled in order to obtain a more crystalline and therefore lower energy structure. By analogy, 
simulated tempering uses a set of discrete values of a temperature parameter T to label and describe 
natter versions of the original posterior distributions. The value T = 1 is reserved for the the 
original posterior distribution. Higher values of T correspond to flatter distributions. In simulated 
tempering, the distribution is 'warmed up' by increasing T. In these natter versions, it is easier for 
the sampler to jump out of local minima and explore the full posterior to find the global minimum. 
Inferences are drawn from the T = 1 sampler. 

As above, let p{H\D,X) be the target posterior distribution we want to sample; by Bayes' 
theorem 

p(H\B,l) oc p{H\l) x p(D\H,I) (Al) 

where we have dropped the normalization factor l/p(D\I). Other posterior distributions at different 
annealing temperatures (3 = 1/T are constructed as 

vr(F|D,X,/3) = p(H\l)p(T>\H,lf (A2) 

= p(H\l)e W (f3 log [p(D\H,l)}) (A3) 

where < /3 < 1. The parameter (3 varies from to 1; (3 = 1 corresponds to the original, target 
distribution, with lower values corresponding to flatter (higher temperature) versions of the target 
distribution. 

In parallel tempering, multiple MCMC chains are run in parallel at tit temperatures {1, Po, f3i, ..., /3 nT } 
for tit > 1- At intervals, proposals are made to swap the parameter states at adjacent but ran- 
domly selected temperatures. For example, at iteration t, suppose that the sampler at /% has a 
parameter H tt i, and /3j + i has a parameter state H ty i+\. These are the candidate parameter states 
for swapping. The swap is accepted with probability 

. f 7 r(g M+1 |D,/3 i ,X)7r(g M |D,/3 i+ i,X) ) 
r mm \7r(^|D,A,X)7r(F tim |D,/3 m ,X)i- l > 

The swap is accepted if XJ\ ~ Uniform[0, 1] < r, that is, if a number U\ drawn from a uniform 
random distribution between zero and 1, is less than or equal to r. If the swap is accepted, then 
the parameter states are swapped: the chain indexed i now has parameter state H^i+i, and the 
chain indexed i + 1 now has parameter state .fft j. This swapping process propagates information 
across the parallel simulations. At higher temperatures, the algorithm can explore very different 
locations in the posterior parameter space. At lower temperatures, the algorithm can improve 
local knowledge of the space around minima. Swapping allows highly probable parameter states to 
propagate down to lower temperatur es where t hey c an be explored locally. The swap itself need 



not be proposed at every iteration. iGregoryl (J2005I ) implements an example parallel tempering 



algorithm by allowing a swap on average once every n s iterations: the swap is only performed if the 
value of U2, drawn from a uniform distribution between zero and 1, is less than or equal to l/n s . 



Each of the MCMC chains uses the Metropolis-Hastings algorithm ([Gregory||2005l ) to explore 
each 7r(.ff~|D,Z, /3). Normal distributions were used as the proposal distributions for Metropolis- 
Hastings algorithm. Widths for each proposal distribution were found after making several shorter 
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exploratory runs of the j3 = 1 chain with an adaptive algorithm that varied the proposal distribution 



with to generate an acceptance ratio in the range 0.16 — > 0.30 (JGelman et al.ll2003l ). For each 
variable 9 in each spectral model, a uniform prior is assumed, that is, p(6) = l/(#i — #o) for 
#o < 9 < #i and p(6) = otherwise. The lower (Qq) and upper (#i) values are constants. The 
limits(6>o) and upper (0i) and the proposal distribution step-size are given in Tables |U 

As is described in the main text, the parallel tempering MCMC algorithm produces marginal 
distributions of E c for the 23 July 2002 flare consistent with expectations. The parallel tempering 
MCMC algorithm described here was used in the analysis of both the 19 January 2005 and 23 July 
2002 flares. 



B. Implementation of the parallel tempering Markov chain Monte Carlo algorithm 

The results described in the paper arise from implementing the parallel tempering algorithm 
described in Section O Five temperatures in the algorithm are used: f3 = 1,0.75,0.5,0.25,0.01. 
Each simulation takes 50,000 samples (five times as many samples as the Monte Carlo approach 
of Section 13.1.30 . The simulation is run ten times with a different starting point chosen uniformly 
randomly in the volume 6q — 5s, #o+5s, where s is the size of the proposal distribution step size. The 
proposal distribution step size is the square root of the diagonal elements of the covariance matrix 
of a least squares fit calculated at 9q. The last half of the samples are considered post burn- in, 
and are retained. Convergence b e tween and within the 10 simulation runs is assessed using the R 



measurement from lGelman et al.l (|2003l ). In all cases , the .R-measurement was below approxim ately 



1.1, which may be taken as indicating convergence ( Gelman et al.1 120031 ; Ivan Dyk et al.ll200ll ). 



C. Normality of the marginal distributions 

The normality of the univariate marginal distributions was assessed using Q-Q (quantile- 
quantile) plots (Figures [TT1 and [T2]) . A Q-Q plot is a graphical method of comparing two different 
distributions, and is constructed as follows. The cumulative distribution function of a random 
variable X is defined as 

F x (x) = P(X < x) (CI) 

that is, the probability that the random variable X takes on a value less than or equal to x. The 
function Fx{x) is monotonically increasing in the range zero to one. The inverse of Fx is called 
the quantile function, Q, and is defined as 



Qx(r) = x if F x {x) 



(C2) 



If Fx is a one-to-one function, the inverse Q is uniquely determined. If the function Fx is not one- 
to-one the inverse Q can be defined as the weighted average of all relevant points. The definition 
of the quantile function applies to random variables or sample distributions. The Q-Q plots shown 
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Fig. 7. — Results from each of the four uncertainty analysis methods (Section [3j) for (a) EM and 
(b) kT from the model fit to the 23 July 2002 flare data. These plots follow the same convention 
as Figure El See Section S] for more detail on these results. 



Table 4. Details of the prior variable ranges and the proposal distribution step-size used in the 

Bayesian/MCMC analysis of the 19 January 2005 and 23 July 2002 flare data. Priors for each 

variable are uniform within the stated ranges. Each proposal distribution is Normal, with width 

as indicated. See Section [2] for more detail on the choice of model, and Appendix [A] for more 
detail on the implementation of the Bayesian/MCMC analysis, the proposal distributions are all 

Normal. 



Flare 


Parameter 


Prior range 


Proposal distribution width 


19 January 2005 


EM 


0.77 -> 6.94 


0.01 




kT 


0.68 -> 6.08 


0.01 




F 


0.01 -+ 1000 


0.0004 




<h 


1.1 -> 20 


0.002 




E c 


6.8 ->■ 290 


0.17 




(h 


3334 -* 33347 


821 




G 2 


1279 ->• 12790 


283 


23 July 2002 


EM 


0.9 -+8.14 


0.004 




kT 


0.5 ^8.0 


0.001 




A 


0.002 -> 0.3 


0.0003 




Si 


1.1 -> 50 


0.014 




E b 


50 -> 32000 


7.5 




<h 


1.1 -» 50 


0.007 




E c 


0.01 -> 50 


0.57 
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Fig. 8. — Results from each of the four uncertainty analysis methods (Section [3]) for (a) A, (b) 5\, 
(c) Eh, (d) 62 and (e) E c , from the model fit to the 23 July 2002 flare data. These plots follow the 
same convention as Figure El See Section 0] for more detail on these results. 
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Fig. 9. — Two dimensional marginal probability density functions for the parameters of the model 
used to fit the spectrum of the 23 July 2002 flare. In contrast to similar distributions plotted in Fig. 
[5] for the 19 January 2005 flare, some distributions are highly asymmetric within the parameter 
ranges found. The number on the upper right of each plot is the Spearman rank correlation 
coefficient for the abscissa versus the ordinate. There are many more moderately and strongly 
(anti-) correlated pairs of parameters for this flare model compared to the 19 January 2005 flare 
model. For some pairs of parameters (for example 5\ versus A and 5^ versus 2£&), the proportion of 
the space taken up by the high probability volume is relatively small, and for others (for example, 
E c versus A), it is relatively large. For the model applied to this flare spectrum, many of the 
resulting probability density functions do not show Normal distribution shapes. This indicates 
that the hypersurface for the model fit to these flare data has a more complicated structure than 
the hypersurface of the model fit to the 19 January 2005 flare. 
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Fig. 10. — Electron spectrum results for the flare- injected electrons arising from the 
Bayesian/MCMC method for the 23 July 2002 flare, (a) Electron spectrum (flux (in units of 
erg keV - s _1 ) multiplied by E 3 ' 38 ) with 68% and 95% credible interval spectra indicated by the 
dashed and dotted lines, respectively. The electron flux spectrum corresponding to Q MAP is indi- 
cated by the solid line, (b) 68% and 95% credible intervals (dashed and dotted lines, respectively) 
relative to the Q MAP electron flux spectrum, (c) Flare injected electron number flux probability 
density function, with 68% and 95% credible intervals indicated, (d) Flare injected electron power 
probability density function, with 68% and 95% credible intervals indicated. In plots (c) and (d) 
the distribution mean/mode is indicated by the solid/dot-dashed vertical line. 
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in Figures [TT] and [12] show a set of open circles and a straight line. A circle is plotted at the point 
where the abscissa and the ordinate are the values of quantile functions for the standard Normal 
distribution N(0, 1) and the marginal distribution, for a given value of probability r. A straight line 
is drawn through the points defined by the quantile functions for the standard Normal distribution 
and a Normal distribution N ( 9i, oj. ) , where 



N S 



N s . 



and 

N S 
_ 

~ N s 



^E{i-^ 



are estimated from the Ns samples of the parameter Oi, 1 < i < N$. The straight line enables an 
assessment of how closely the marginal distribution follows a Normal distribution, and where any 
deviations occur. The quantile function for the standard Normal distribution function is called the 
probit function and is defined as 

probit(r) = y/2 erf" 1 (2r - l),r € (0, 1) (C3) 

where erf - (x) is the inverse error function. The probit function gives the value of a iV(0, 1) random 
variable associated with specified cumulative probability r, for example: 

probit(0.025) ~ -1.96 ~ -probit(0.975). (C4) 

Therefore, and conveniently, the abscissa in the Q-Q plots can be understood as multiples of the 
standard deviation away from the mean. The Q-Q plots were implemented usi ng the 'R' statistical 



comp uting environment, available from the R Project for Statistical Computing (JR Development Core Team 



2011 
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Fig. 11. — Q-Q plots for the Bayesian/MCMC samples of the 19 January 2005 model spectrum 
parameter values. All parameters are approximately Normally distributed in the range —2, +2 
quantiles about the estimated mean. The tails of the distributions show deviations away from a 
true Normal distribution. Curvature of the sample distribution at negative quantiles indicates that 
the tail is thinner than that expected from the sample Normal distribution N ( $i, o~q. ) . Similarly, 
curvature of the sample distribution at large positive quantities indicates that the tail is fatter than 
that expected from the sample Normal distribution N (#j,(70 
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Fig. 12. — Q-Q plots for the Bayesian/MCMC samples of the 23 July 2002 model spectrum param- 
eter values. The two thermal parameters of the model EM and kT appear to be approximately 
Normally distributed; the remaining non-thermal parameters (A, S±, E^, 62 and E c ) are clearly not 
Normally distributed. 
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